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Abstract 

An important task of uncertainty quantification is to identify the probability of 
undesired events, in particular, system failures, caused by various sources of uncer¬ 
tainties. In this work we consider the construction of Gaussian process surrogates 
for failure detection and failure probability estimation. In particular, we consider 
the situation that the underlying computer models are extremely expensive, and in 
this setting, determining the sampling points in the state space is of essential impor¬ 
tance. We formulate the problem as an optimal experimental design for Bayesian 
inferences of the limit state (i.e., the failure boundary) and propose an efficient 
numerical scheme to solve the resulting optimization problem. In particular, the 
proposed limit-state inference method is capable of determining multiple sampling 
points at a time, and thus it is well suited for problems where multiple computer 
simulations can be performed in parallel. The accuracy and performance of the 
proposed method is demonstrated by both academic and practical examples. 
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1 Introduction 


Real-life engineering systems are unavoidably subject to various uncertainties 
such as material properties, geometric parameters, boundary conditions and 
applied loadings. These uncertainties may cause undesired events, in particu¬ 
lar, system failures or malfunctions, to occur. Accurate identihcation of failure 
region and evaluation of failure probability of a given system is an essential 
task in many helds of engineering such as risk management, structural design, 
reliability-based optimization, etc. 

Conventionally the failure probability is often computed by constructing linear 
or quadratic expansions of the system model around the so-called most proba¬ 
ble point, known as the hrst/second-order reliability method (FORM/SORM); 
see e.g., [T5II3U] and the references therein. It is well known that FORM/SORM 
may fail for systems with nonlinearity or multiple failure regions. The Monte 
Carlo (MC) simulation, which estimates the failure probability by repeatedly 
simulating the underlying system, is another popular method for solving such 
problems. The MC method makes no approximation to the underlying com¬ 
puter models and thus can be applied to any systems. On the other hand, 
the MC method is notorious for its slow convergence, and thus can become 
prohibitively expensive when the underlying computer model is computation¬ 
ally intensive and/or the system failures are rare and each sample requires to 
a full-scale numerical simulation of the system. To reduce the computational 
effort, one can construct an computationally inexpensive approximation to the 
true model, and then evaluate the approximate model in the MC simulations. 
Such approximate models are also known as response surfaces, surrogates, 
metamodels, and emulators, etc. These methods are referred to as the re¬ 
sponse surface (RS) methods |20lll0|l37ll21ll22ll36j in this work. The response 
surface can often provide a reliable estimate of the failure probability, at a 
much lower computational cost than direct MC simulations. 

In this work we are focused on a specihc kind of RS, the Gaussian processes 
(GP) surrogates, also known as kriging in many helds of applications. The 
GP surrogates have been widely used in machine learning |31], geostatistics 
[31], engineering optimizations [ID], and most recently, uncertainty quantihca- 
tions m- A number of GP-based methods have been also been successfully 
implemented for failure probability estimation [28E] . In this work we consider 
the situation where the underlying computer models are extremely expensive 
and one can only afford a very limited number of simulations. In this setting, 
choosing the sampling points (i.e. the parameter values with which the simu¬ 
lation is performed) in the state space is of essential importance. Determining 
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the sampling points for GP can be cast as to optimally design computer exper¬ 
iments and considerable efforts |2S] have been devoted to it. Those methods 
aim to construct a surrogate that can accurately approximate the target func¬ 
tion in the whole parameter space. As will be explained later, in the failure 
probability estimation or failure detection problems, only the sign of the tar¬ 
get function is used. Thus by requiring surrogates to be globally accurate, the 
methods may allocate considerable computational efforts to the regions not of 
interest, and use much more model simulations than necessary. 

Several methods have developed to determine sampling points for the failure 
probability estimation purposes. Most of these methods consist of sequentially 
hnding the “best point” as a result of a heuristic balance between predicted 
closeness to the limit state, and high prediction uncertainty, e.g. [T71I7] . Such 
methods are shown to be effective in many applications, while a major lim¬ 
itation is their point-wise sequential nature, which makes it unsuitable for 
problems in which multiple computer simulations can be performed parallelly. 
An exception is the stepwise uncertainty reduction (SUR) method developed 
in [SHE], in which the authors proposed an optimal experimental design frame¬ 
work which determines multiple sampling points by minimizing the average 
variance of the failure probability. It should be noted that the design criteria 
in the SUR method is particularly developed for the goal of estimating the 
failure probability only. In practice, one is often not only interested in esti¬ 
mating the failure probability, but also identifying the events that can cause 
failures; the latter demands a design criteria for the goal of detecting the limit 
state, i.e., the boundaries of the failure domain. In this work, we recast the 
surrogate construction as a Bayesian inference to identify the limit state, and 
based on that, we formulate an information-theoretic optimal experimental 
design, which uses the relative entropy from the prior to the posterior as the 
design criteria, to determine the sampling points. We also present an efficient 
numerical scheme for solving the resulting optimal design problem, modified 
from the simulation-based method developed in |2lj. We compare the perfor¬ 
mance of the proposed limit-state inference (LSI) method with that of the 
SUR by numerical examples. 

We note that another line of research in failure probability estimation is to 
develop more efficient sampling schemes, such as the subset simulations m, 
imoortance samolinEr |321I19]. seauential Monte Carlo nn, the cross-entropy 
method [38IIT3] . etc. For practical engineering systems, computer simulations 
can be extremely time consuming. In many cases, one can only afford very lim¬ 
ited number of simulations — nothing beyond a few hundreds. In this case, 
even the most effective sampling method is not applicable. To this end, surro¬ 
gates are needed even in those advanced sampling schemes and in particular 
the LSI method can be easily integrated into the aforementioned sampling 
schemes, resulting in more efficient estimation schemes. Examples of combin¬ 
ing surrogates and efficient sampling schemes include PTIEHIITIITB] . 


3 









The rest of this paper is organized as following. We hrst review the prelim¬ 
inaries of onr work in Section including the mathematical formulation of 
failure probability computation and the GP surrogates. Our Bayesian experi¬ 
mental design framework and its numerical implementations are presented in 
Section]^ Numerical examples are presented in Section]^ to demonstrate the 
effectiveness of the proposed method, and hnally Section offers some closing 
remarks. 


2 Problem formulation 


2.1 Failure detection and failure probability estimation 


Here we describe the failure probability estimation problem in a general set¬ 
ting. We consider a probabilistic model where x is a d-dimensional random 
variable that represents the uncertainty in model and the system failure is 
often dehned using a real-valued function g{-) : —)■ i?, which is known as 

the limit state function or the performance function. Specihcally, the event of 
failure is dehned as ^'(x) < 0 and as a result the failure probability is 


P = P(fif(x) < 0) = / 4(x)p(x)(ix = f 4(x)p(x)cix, 


'{xeR‘^|g(x)<0} 

where /^(x) is an indicator function: 


fl iff7(x)<0 

[O iff7(x)> 


and p(x) is the probability density function (PDF) of x. In what follows we 
shall omit the integration domain when it is simply P'^. This is a general 
dehnition for failure probability, which is used widely in many disciplines in¬ 
volving reliability analysis and risk management. P can be computed with the 
standard Monte Carlo estimation: 

( 1 ) 

i=i 

where samples xi, ...,x„ are drawn from distribution p(x). 


In practice, many engineering systems require high reliability, namely the fail¬ 
ure probability P <C 1. In this case, MC requires a rather large number of 
samples to produce a reliable estimate of the failure probability. For exam¬ 
ple, for P 10“^, MC simulation requires 10® samples to obtain an estimate 
with 10% coefficient of variation. On the other hand, in almost all practical 
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cases, the limit state function g{x.) does not admit analytical expression and 
has to be evaluated through expensive computer simulations, which renders 
the crucial MC estimation of the failure probability prohibitive. To reduce the 
number of full-scale computer simulations, one can construct a computation¬ 
ally inexpensive surrogate G(x) to replace the real function g(x) in the MC 
estimation. In this work we choose the Gaussian Process surrogates and we 
provide a description of GP in the next section. 


2.3 Gaussian process surrogates 


The GP method constructs the approximation of g{x) in a nonparameteric 
Bayesian regression framework [33IH1] . Specihcally the target function g{x) is 
cast as 

g{x) = /x(x) + e(x) (2) 

where /i(x) is a real-valued function and e(x) is a zero-mean Gaussian process 
whose covariance is specihed by a kernel K{x,x'), namely, 

COV[e(x),e(x')] = K{x,x). 

The kernel K (x, x') is positive semidehnite and bounded. Suppose that N 
computer simulations of the function g{x) are performed at parameter values 
X* := [x^,. .. X*], yielding function evaluations y* := [yl,.. .y^], where 

y* = 9{^i) for i = 

Suppose we want to predict the function values at points D := [xi,.. .x^], 
i.e., y = [ 2 / 1 ,... 2/m] where yi = g(xi). The posterior distribution of y is 

y I D,X*,y*~iA^(u,C), (3a) 

where the posterior mean u = [ui,... Um) is a m-dimensional vector with each 
element to be 


Uj = y{xj) + rjR ^(y* — pt) for j = l...m, (3b) 

where pL = {pi, ..., pn) and pi = pix*) for i = 1,... ,n. The posterior covari¬ 
ance matrix C is given by 

(C )jj' = K{xj,xp) - rj R-^ rp. (3c) 

In the equations above, rj represents a n-dimensional vector whose Tth com¬ 
ponent is K(x*,Xj) and the matrix R is given by (R)*,*' = K(x*,x*,) for 
i, i' = 1,..., n. 
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3 The experimental design framework 


A key question in constructing a GP surrogate is to determine the locations D 
where the true limit state function g{x.) is evaluated, which is often known as 
a design of computer experiments. In this section, we show that determining 
the sampling locations can be translated into a Bayesian experimental design 
problem whose goal is to hnd the locations of limit state with a given prior 
distribution. 

3.1 Bayesian inference experiments 

We formulate a Bayesian inference problem based on the following argument. 
In failure probability estimation, the limit state function g is only used in 
the indicator function /g(x) and so one is really interested in the sign of ^'(x) 
rather than the precise value of it. To this end, the essential task in con¬ 
structing surrogates for the failure probability estimation is to learn about 
the boundary of the failure domain. Let z represents the boundary of the fail¬ 
ure domain, i.e., the solutions of g{x.) = 0. We want to identify the boundary 
z with Bayesian inference. In this setting, the data is obtained by making 
observations, i.e., evaluate g{x.) at locations D = (xi,X2,... ,x„), resulting in 
data y = {yi,... ,yn) with each yi = 5 '(xj). The likelihood function p(y|z,D) 
in this case is given by Eqs. 0 with X* = [z] and y* = [0]. We also need 
to choose a prior distribution for z, and without additional information, it is 
reasonable to assume that the prior is simply p(z). In this setting the posterior 
distribution of z with data y can be computed with the Bayes’ theorem: 



where p(y|D) is the evidence. 

3.2 Optimal design criteria 

As is mentioned earlier, determining the locations D can be formulated as a 
design of the Bayesian inference experiments. Following the decision-theoretic 
optimal experimental design framework, an objective for experimental design 
can be generally formed: 



J J ui'D,y,z)p{z\y,'D)p{y\B)dzdy, (4) 


6 



where m(D, y, z) is a utility function, U (D) is the expected utility. The utility 
function u should reflect the usefulness of an experiment at conditions D, 
i.e., we will get a particular value outcome y at condition D by inputting a 
particular value of the parameters z. Since we do not know the exact value of 
z and y in advance, we take the expectation of u over the joint distribution of 
z and y, resulting in the expected utility ?7(D). The optimal choice of D then 
can be obtained by maximizing the expected utility of the design space S: 

D* = argmaxt/(D). (5) 


A popular choice for the utility function is the relative entropy, also known as 
the Kullback-Leibler divergence (KLD), between the prior and the posterior 
distributions. For two distributions pyi(z) and ps(z), the KLD from pA to pb 
is dehned as 


Dkl{pa il Pb) = f Pa{z) ln[ ^^|^| ]dz = E^pn (6) 

J Pb(z) Pb(z) 

where we dehne OlnO = 0. This quantity is non-negative, non-symmetric, and 
reflects the difference in information carried by the two distributions. When 
KLD is used in our inference problem, the utility function becomes 


M(D,y,z) = ®;^i(p(-|y,D) || p(-)) = / p(z|y,D)ln 


p(z|y,D) 

p(z) 


dz. (7) 


The utility function m(D, y, z) in Eq. ([^ can be understood as the information 
gain by performing experiments under conditions D, and larger value of u 
implies that the experiment is more informative for parameter inference. Note 
that the utility function in Eq. Q is independent of z and as such the expected 
utility U (D) is reduced to 


U{D) = J ®xi(p(-|y,D) II p{-))p{y\'D) dy 

(z|y,D)ln 


p(z|y,D) 

P(z) 


dzp(y|D)dy, (8) 


where z in Eq. ([^ is replaced by z for simplicity. Next we discuss how to 
numerically solve the optimization problem Eq. ([^, with t7(D) is given by 
Eq. 


3.3 Numerical implementation 


Following the recommendation of [21], we use the simultaneous perturbation 
stochastic approximation (SPSA) method, to solve the optimization prob- 
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lem ([^. SPSA is a derivative-free stochastic optimization method that was 
first proposed by Spall nmg, and we provide the detailed algorithm of SPSA 
in Appendix A. Note here that, since it is a derivative-free method, the algo¬ 
rithm only uses the function value of 17(D). 


Next we discuss the evaluation of 17(D). To start, we re-write Eq. (|^ as 

p(z|y,D) 


17(D) = 


z|y,D)ln 


p(z) 


clzp(y|D) dy 


= / In |C| p(z) dz- J J ln[p(y |D)]p(y |D) dy + Z 


= -2E[ln|C|]-fe[p(y|D)]-f Z 


(9) 


where ]HI(-) is the notation for entropy and Z is a constant independent of D. 
The detailed derivation of Eq. (|^ is given in Appendix B. Note that Eq. ([^ 
typically has no closed-form expression and has to be evaluated with MC 
simulations. Draw m pairs of samples {(yi, Zi),..., (y^, z^)} from p(y, z|D), 
and the MC estimator of E[ln |C|] is 


-1 m 

(7i = -^|C(z„D)|. 


( 10 ) 


Recall that in [2l], the entropy term ]HI[p(y|D)] is computed using a nested 
MC. For efficiency’s sake, we use the resubstitution method developed in pQ to 
estimate the entropy. The basic idea of the method is rather straightforward: 
given a set of samples {yi,..., y^} of p(y|D), one hrst computes an estimator 
of the density p(y|D), say p(y), with certain density estimation approach, and 
then estimates the entropy with. 
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1 

m 


Ep(y*)- 

i=l 


( 11 ) 


Theoretical properties of the resubstitution method are analyzed in [2^23] and 
other entropy estimation methods can be found in [6] . In the original work [1] , 
the distribution p is obtained with kernel density estimation, which can become 
very costly when the dimensionality of y gets high, and to address the issue, 
we use Gaussian mixture based density estimation method 1311 . Without loss 
of generality we can simply set the constant Z = 0 and thus an estimate of U 
is simply 

17 = 17i + 172, (12) 

which will be used to provide function values in the SPSA algorithm. Finally 
we reinstate that the numerical implementation in this work differs from a 
more general implementation outlined in [23] in the following two aspects: 
hrst, thanks to the Gaussianity of p(y|D) , we can analytically integrate the 
variable y in the hrst integral on the right hand side of Eq. 0 ; second, we 







use Gaussian mixture density estimation method to approximate p(y|D) to 
avoid the nested MC integration. We note that there are other methods such 
as [in] that can be used to approximate the nested integral. 


3.4 Sequential optimal design 


Until this point we have presented our method under the assumption that the 
sampling points are determined all in advance of performing computer simula¬ 
tions, which is often referred to as an open-loop design. In many applications, 
a more practical strategy is to choose the sampling points in a sequential 
fashion: determine a set of sampling points, perform simulations, determine 
another set of points based on the previous results, and so forth. A sequen¬ 
tial (close-loop) design can be readily derived from the open-loop version. For 
computational efficiency, we adopt the “greedy” approach for the sequential 
design, while noting that this approach can only obtain sub-optimal solutions 
in general. Simply speaking, the greedy sequential design iterates as follows 
until a prescribed stopping criteria is met: 

(1) determine n sampling points with an open-loop design; 

(2) perform simulations at the chosen sampling points; 

(3) use the simulation results to update p(z) and p(y|z, D); 

(4) return to step (1). 

A key ingredient in the greedy procedure is that the prior p(z) and the likeli¬ 
hood p(z|z, D) are updated based on the simulation results in each stage. To 
be specihc, the prior distribution at step A; -|- 1 is the posterior distribution of 
z at step k\ 

Pfc+i(z) =pfc(z|y,D). 

The update of the likelihood is a bit more complicated. Namely, let repre¬ 
sent all the sampling points at which the function g{-) has been evaluated in 
the hrst k stages and y^ be the associated function values, and the likelihood 
function pfc+i(y|z,D) is once again given by Eqs. ([^ while 

X* = [Xfc,z] and y* = [yfc,0]. 

Note that we do not have analytical expression for pfc_|_i(z) or equivalently 
Pfc(z|y,D); however, our implementation only requires samples drawn from 
Pfc+i(z) which are available from the previous iteration. Another issue that 
should be noted is that, every time a new observation (z,0) is added to X*, 
the posterior covariance needs to be updated. Recomputing the covariance 
matrix for each sample can be rather time consuming, and so here we use the 
covariance matrix update method which computes the new covariance matrix 
based on the previous one. The detailed update formulas are given in [T511T5] 
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Fig. 1. The rescaled Branin function. The solid line is the limit state g{x.) = 0. 

and we shall not repeat it here. Finally we also need to specify a stopping cri¬ 
teria for the algorithm and following the previous works we stop the iterations 
when the difference the estimated probabilities in two consecutive iterations 
is smaller than a threshold value. 


4 Numerical examples 


4-1 Branin function 


To compare the performance of our method with SUR, we hrst test both 
methods on the rescaled Branin function 

g{xi,X2) = 80 - [(15x2 - A(15a;i - 5)^ + ^ (15xi - 5) - 6)^ 

tt 

+ 10(l-^)cos(15xi-5) + 10]. (13) 
Stt 

which is plotted in Fig. Here we assume that Xi and X 2 are independent and 
each follows a uniform distribution in [0,1]. The function is included in the R 
package Kriginv na as an example for the SUR method and the numerical 
results suggest that SUR performs well for this numerical example. 

We hrst perform a standard MC with 10,000 samples, resulting in a failure 
probability estimate of 0.256, which will be used as the “true” failure prob¬ 
ability to validate our results. In the experimental design, we hrst choose 4 
points with the Latin hypercube method as the initial design points for both 
methods. We then sequentially choose 9 points with both methods, where one 
point is determined in each iteration. In the numerical experiment, we choose 
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Fig. 2. Left: the probability of detection error plotted against the number of de¬ 
sign points. Right: the error in the failure probability estimates plotted against the 
number of design points. 

to use a squared exponential covariance: 


K(x,x') 


a exp( 



2 



(14) 


where we choose a = 0.1 and /? = 1 in this example. If desired, the parameter 
values in Eq. (14) can also be obtained with maximum marginal likelihood 
method. The prior mean is determined with a quadratic regression. In all the 
three examples, the algorithm is terminated after 50 iterations in the SPSA 
method, and in each iteration, 5000 samples are used to evaluate the value of 
f/(D). We draw 10000 samples from the uniform distribution, and each time a 
new design point is added, we use the resulting GP surrogates to detect which 
points are in the failure region and which are not. Using the detection results 
we compute the failure probability. To compare the performance of the two 
methods, we compute two quantities as the indicators of performance. The 
Erst is the error between the failure probability estimate by each method and 
the true failure probability which is estimated by standard MC. The second 
is the probability that a point is mis-identihed: 


P[x G {x|^(x) < 0,^(x) > 0} U {x|^(x) > 0,^(x) < 0}], 

where g(x) represents the surrogate. In Fig. we plot both indicators as a 
function of the number of design points for our method and the SUR. The 
results of both indicators suggest that our LSI method seems to have a better 
performance than SUR in this example. 


4-3 Four branch system 


Our second example is the so-called four branch system, which is a popular 
benchmark test case in reliability analysis. In this example the limit state 
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Fig. 3. The limit state function of the four-branch model. The solid lines in both 
figures are the limit state ( 7 (x) = 0. 



Fig. 4. Left: the probability of detection error plotted against the number of de¬ 
sign points. Right: the error in the failure probability estimates plotted against the 
number of design points. 

function reads 


3 -b 0.1(xi - X2)‘^ - (xi -b X2)/\/2 


g{xi, X 2 ) = min < 


3 -b 0.1(a:i - X2)^ + (xi -b X2)/V^ 
{xi - X 2 ) + 7/a/2 




^ {X2 - Xi) - 7/y/2 


which is shown in Fig. The input random variable Xi and X 2 are assumed to 
be independent and follow standard normal distribution. We hrst compute the 
failure probability with a standard MC estimation of 10^ samples, resulting 
an estimate of 2.34 x 10“^, i.e., 234 samples fall in the failure region. 


We use the squared exponential covariance function (14) in this example, but 
with a = 0.2 and {3 = 0.8. The prior mean of the GP surrogate is computed 
with a quadratic regression. We perform the sequential design described in 
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Fig. 5. Top: the left figure shows the initial design points and the right one shows the 
design points determined in the first iteration. Bottom: the left and the right figures 
are the design points determined in the 7th and the 13th iteration respectively. 


Section with 4 sampling points determined in each iteration. The algorithm 
terminates in 13 iterations, resulting in totally 62 design points. We plot the 
errors in failure probability estimation as a function of the number of iterations 
in Fig. 1^ (left), and the mis-detection probability, also as a function of the 
number of iterations in Fig. (right). In Figs. we plot the initial design 
points and the design points found in the hrst, the 7th, and the last iterations. 
In Fig.|^(left), we show all the 62 design points including both the initial points 
and those found by our method. Also shown in the hgure is the surrogate 
constructed with the obtained design points. We can see that our method 
allocates more points near the boundary of the failure region than random 
sampling points. In Fig. (right), we show the samples that are incorrectly 
identihed (crosses) and those are correctly identihed as failures (dots). 


We note that, the failure probability estimate obtained in the hnal iteration is 
2.31 X 10“^, while the estimate of standard MC is 2.34 x 10“^. Also, compared 
to the MC simulation with true model, we hnd that totally 19 samples are 
incorrectly classihed: 11 safe samples are identihed as failures, and 8 failure 
samples are identihed as safe ones. The numerical results indicate that our 
method can obtain reliable estimates of the failure probability with a rather 
small number of sampling points. 
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Fig. 6. Left: all the design points (asterisks), the limit stage (solid line) and the GP 
surrogate (dashed line). Right: the dots are the points that are corrected identified 
as failures and the crosses are the mis-identified samples. 

4-3 Clamped beam dynamics 


Our last example is a practical problem which concerns the dynamic behavior 
of a beam clamped at both ends and a uniform pressure load is suddenly 
applied as shown in Fig.[^ We are interested in the dynamics of the deflections 
at the center of the beam caused by the pressure. 


The Lagrangian formulation of the equilibrium equations without body forces 
can be written 


d^u 


V • (F ■ S) 


0 , 


where p is the material density, u is the displacement vector, F is the defor¬ 
mation gradient, 


F = 


duj 

dxk 


+ di 


ik 


and S is the second Piola-Kirchoff stress tensor. For simplicity, we assume 
linear elastic constitutive relations and isotropic material. As a result, the 
constitutive relations may be written in matrix form: 


✓ N 

^11 


1 

(M 

0 

0 

1_ 


f \ 

Ell 

S22 

> = 

C12 C22 

< 

E22 

Sl 2 

K J 


1 

bO 

to 

1_ 


Ei 2 


where 


and 


1 f ^ duk duk \ 

2 \dxj dxi dxi dxj J ’ 


Cii — — 


22 



Cl2 — 


Eu 

1 - 1 / 2 ’ 


G 


12 — 


E 

2(1 + z/)‘ 
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parameter 

L 

h 

P 

E 

V 

(5 

mean 

5 

0.1 

7.77"^ 

3 X 10'^ 

0.3 

200 

variance 

6.25 X 10-2 

2.5 X 10“® 

1.51 X 10-9 

3.6 X 10^® 

2.25 X lO-'^ 

40 


Table 1 

The mean and variance of the random parameters in the clamped beam problem. 
Here E is Young’s modulus and u is Poisson’s ratio. The initial conditions are 

m(0,x) = 0, ^(0,x)=0. 

Readers who are interested in more details about the Lagrangian formulation 
for nonlinear elasticity can consult, for example, [30] . 

We assume in this examples that the beam length L, the height h, the material 
density p, the Young’s module E, the Poisson ratio z/, and the applied load S 
are random. All the random parameters follow a normal distribution and are 
independent to each other. The means and the variances of the parameters 
are summarized in Table 1. To demonstrate the statistical variations of the 
beam dynamics, we randomly generate 10 parameter samples and plot the 
10 resulting beam dynamics in Fig. We dehne the failure as the maximum 
deflection at the center of the beam being larger than a threshold value Umax- 
To be specihc, we take Umax = 0.23. We hrst run standard MC with 10® samples 
to compute the failure probability, resulting in an estimate of 3.35 x 10“^. 

In the GP method, we use 64 initial sampling points drawn by the hyper 
Latin cube approach. Our algorithm determines 4 points in each iteration and 
it is terminated after 25 iterations, resulting in totally 164 sampling points. 
As before, we plot the errors in failure probability estimation and the mis- 
detection probability, both against the number of iterations, in Figs. In 
both plots we can see that the accuracy of the estimations improve rapidly 
as the number of sampling points increases, indicating that our methods can 
effectively identify good sampling points for this practical problem. Our failure 
probability estimate is 3.35 x 10“^ while 12 points in the failure region is mis- 
identihed as safe ones and 10 safe ones are mis-identihed as failures. 


5 Conclusions 


In conclusion, we have presented a GP-based failure-detection method to 
construct GP surrogate for failure probability estimations. In particular, the 
method recasts the failure detection as inferring a contour g{x.) = 0 with 
Bayesian methods, and then determines the optimal sampling locations for 
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Fig. 7. Schematic illustration of the clamped beam. 



Fig. 8. Dynamics of the deflection at the beam center for ten randomly generated 
samples. 




Fig. 9. Left; the probability of detection error plotted against the number of de¬ 
sign points. Right: the error in the failure probability estimates plotted against the 
number of design points. 
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the inference problem. An efficient numerical implementation of the result¬ 
ing design problem, based on the SPSA algorithm, is also presented. With 
numerical examples, we demonstrate that the proposed LSI method can ef¬ 
fectively and efficiently determine sampling points for failure detection and 
failure probability estimation. 

There are some improvements and extensions of the method that we plan to in¬ 
vestigate in the future. First, just like any GP-based method, one must choose 
the covariance kernels, the prior mean functions, and various hyperparame¬ 
ters. We hope to develop effective methods to determine those functions and 
parameters specially for the failure probability estimation problems. Secondly, 
as we have mentioned, the method shares a lot of common features of the SUR 
approach; comprehensive comparison of the performance, and detailed analysis 
of the advantages and limitations of the two methods, are a very interesting 
problem to us. Other possible improvements include to replace the greedy 
method with a dynamical programming approach in the sequential design, 
and to use surrogates that combine both GP and polynomial chaos. Finally, 
we note that, in addition to failure probability estimations, the LSI method 
can be applied to many other problems as well. In particular, we plan to use 
the method to approximate the feasible regions in constrained optimization 
problems |35j . 
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A The SPSA method 


Here we briefly introduce the SPSA method, in the context of our specihc 
applications. In each step, the method only uses two random perturbations 
to estimate the gradient regardless of the problem’s dimension, which makes 
it particularly attractive for high dimensional problems. Specihcally, in step 
j, one hrst draw a dimensional random vector Aj = [Sj^i, where 

Ud = dim(D), from a prescribed distribution that is symmetric and of hnite 
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inverse moments. The algorithm then updates the solution using the following 
equations: 


Dj+i — Dj 

C/(Dj+c,A,)-(7(D,-c,A.' 



(A.lb) 


(A. la) 


where f/(-) is computed with Eq. (12) and 



with A,a,c, and 7 being algorithm parameters. Following the recommendation 
of [cite], we choose Aj ~ Bernoulli(0.5) and A = 100, a = 0.602, c = 1, and 
7 = 0.101 in this work. 


B Derivation of Equation (3.6) 

From Eq. (|^, we have 



Recall that p(y|z,D) is multivariate normal distribution: 



Because u and C only depend on z and D ( /i = ^{z,d),C = C{z,d) ), we 
change of variable: 


s = C-'/yy - u) 
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1 


ln(27r) 


1 

2 


In |C|(iz 


s's 


= —- / In \C\p{z)dz — I ( -s's 


j^^exp{-ys)p{z)dzds 


1 \ Tl 

7 —r—77 exp(— s's)ds -ln(27r) 

(27r)’^/2 2 ^ 2 ^ ^ 


(B.2) 


Note that the second integral on the right hand side of Eq. (B.2) actually does 
not depend on D and so we can dehne a constant Z such that 


Z = - 

It follows immediately that 


-s s 
2 


1 / 1 / \ , ’ 2 , 

7 —exp(—s sjds-ln(27r) 

(2vr)V2 2 ^ 2 ^ ^ 


ln(p(y|z,D))p(z|z,D)p(z)cizdy = --E;,[ln|C|] + Z, 


which in turns yields Eq. (§. 
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